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ABSTRACT 


Landing  aircraft  on  board  carriers  is  a  most  delicate  phase  of  flight  operations  at 
sea.  The  ability  to  predict  the  aircraft  carrier's  motion  over  an  interval  of  several  sec¬ 
onds  within  reasonable  error  bounds  may  allow  improvement  in  touchdown  dispersion 
and  a  more  certain  value  for  a  ramp  clearance  due  to  a  smoother  aircraft  trajectory. 
Also,  improved  information  to  the  Landing  Signal  Officer  should  decrease  the  number 
of  waveofTs. 

This  work  indicates  and  shows  graphically  that,  based  on  the  data  for  pitch,  heave 
and  roll  measured  for  various  ships  and  sea  conditons,  the  motion  can  be  predicted  well. 
The  predictor  was  designed  on  the  basis  of  Kalman's  optimum  filtering  theory  for  the 
discrete  time  case,  adapted  for  real-time  digital  computer  operation. 
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I.  INTRODUCTION 


The  landing  phase  of  VTOL,  STOL  and  HELO  aircraft  aboard  a  ship  without  a 
ship's  Inertial  Navigational  System  represents  a  complex  operation  and  a  demanding 
task.  The  last  10  to  15  seconds  before  the  aircraft  touchdown  involves  terminal  guidance 
and  control  problems,  where  not  only  the  aircraft  is  disturbed  by  several  kinds  of 
stochastic  disturbances  (wind),  but  also  the  touchdown  point  (on  the  ship)  is  being 
moved  randomly.  Despite  the  wind  disturbances  and  the  final  point  (target)  random 
motion,  the  landing  accuracy  specified  for  ship  operations  is  very  high,  i.e.,  a  few  tens 
of  feet  longtitudinal  landing  dispersion.  Such  a  terminal  point  problem  is  made  tractable 
in  a  most  natural  way  by  assuming  that  the  ship's  position  can  be  predicted  for  several 
seconds  ahead  so  that  the  aircraft  is  guided  toward  the  future  position  of  the  touchdown 
point.  The  scope  of  this  study  was  to  establish  to  what  extent  a  stochastic  process,  like 
the  ship's  motion,  is  predictable  over  moderate  periods  of  time. 

Graphical  results  obtained  throughout  this  estimator's  feasibility  study  concerning 
the  relationship  between  the  ship's  states:  heave,  pitch  and  roll  versus  the  ship's  esti¬ 
mated  states;  heave,  pitch  and  roll;  and  the  influence  of  noise  measurement  are  pre¬ 
sented.  Digital  simulations  show  that  the  prediction  accuracy  was  very  accurate  for 
different  swell  wave  height  and  period.  The  feasibility  of  estimating  the  ship's  motion 
of  acceptable  bounds  of  error  can  also  lead  to  improvement  of  the  Landing  Signal  Offi¬ 
cer's  (LSO)  decision  policy  for  waveoffs. 

This  study  was  based  on  the  use  of  the  ship's  motion  equations  and  the  simulated 
measuring  instrumentations  existing  on  board  ship  in  order  to  get  the  predicted  motions. 
Using  this  information,  a  predictor  based  on  Kalman's  theory  of  optimum  estimation 
was  designed. 

Several  circumstances  contributed  to  the  success  of  this  approach.  The  size  and 
mass  of  the  ship  significantly  filter  the  motion  of  the  sea.  A  complete  landing  operation 
is  short  enough  that  the  stochastic  processes  are  reasonably  taken  to  be  stationary. 
Finally,  the  prediction  interval  is  only  a  small  fraction  of  the  time  it  takes  each  aircraft 
for  final  approach. 

This  work  is  divided  into  three  parts:  l)  the  derivation  of  the  mathematical  model 
of  the  ship's  motion,  2)  the  rationale  in  the  implementation  of  Kalman  filter  and  pre¬ 
dictor  equations,  and  3)  discussion  of  the  simulated  computer  results  obtained.  Since 


1 


we  are  interested  only  in  the  most  critical  aspects  of  the  landing  operation,  namely  the 
characteristics  of  the  longtitudinal  channel,  we  merely  investigate  in  the  sequel  the  pre¬ 
dictability  of  the  heave,  pitch  and  roll  motion  of  the  ship  where  yaw  is  omitted. 
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II.  PROBLEM  STATEMENT 


A.  GENERAL 

The  ship's  attitude  estimation  used  in  this  thesis  involves  calculation  of  heave,  pitch 
and  roll  of  the  ship  moving  in  the  general  direction  towards  a  swell  while  recovering  an 
aircraft.  The  motion  of  the  ship  is  given  in  z  ,  0  and  4>  coordinates  as  shown  in  Figure 

1. 

This  problem  will  be  developed  using  state  space  methods.  Given  the  heave,  pitch 
and  roll  motion  (the  measurements)  received  through  the  ship's  sensor  system,  we  are 
interested  in  estimating  the  heave,  pitch  and  roll  of  the  ship.  The  state  variables  for  this 
plant  are  zh  ,  i* ,  0,  ,  0t ,  and  <£,  . 

B.  SHIP  MOTION 

A  ship  moving  on  the  surface  of  the  sea  is  almost  always  in  oscillatory  motion.  The 
different  kinds  of  oscillatory  motions  that  a  ship  experiences  can  be  described  with  the 
help  of  Figure  1,  which  shows  the  six  kinds  of  motion,  three  linear  and  three  rotational 
about  the  principal  axes.  Accordingly, 

The  six  kinds  of  ship's  motions  are: 

1. a  =  surging  -  motion  backwards  in  the  direction  of  the  ship  travel. 

2.  b  =  swaying  --  athwartship  motion  of  the  ship. 

3.  c  =  heaving  —  motion  vertically  up  and  down. 

4.  d  =  rolling  -  angular  motion  about  the  longitudinal  axis.  When  the  ship  rolls,  it 
lists  alternately  from  starboard  to  port  and  then  back  to  starboard. 

5.  e  =  pitching  --  angular  motion  about  the  transverse  axis.  When  a  ship  pitches,  it 
trims  alternately  by  the  bow  and  by  the  stern. 

6.  /  =  yawing  -  angular  motion  about  the  vertical  axis. 
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jr-axis  is  longitudinal  axis 
j'-axis  is  transverse  axis 
z-axis  is  vertical  axis 


Only  three  kinds  of  motion,  namely  heaving,  pitching  and  rolling  (  z  ,  6  ,  <f>  )  are 
purely  oscillatory  motions,  since  these  motions  act  under  a  restoring  force  or  moment 
when  the  ship  is  disturbed  from  its  equilibrium  position.  These  ship  motions  are  the  only 
motions  being  considered  in  this  thesis. 

Although  in  reality  a  ship  experiences  all  six  kinds  of  motion  simultaneously,  only 
one  motion  of  the  three  angular  motions  will  be  treated  at  a  time  in  the  following 
sections.  However,  one  must  bear  in  mind  that  any  one  kind  of  motion  is  not  inde¬ 
pendent  of  the  others;  in  the  first  approximation,  coupling  between  the  motions  is  neg¬ 
lected  in  order  to  simplify  the  problem. 

C.  SYSTEM  MODEL 

The  system  to  be  modeled  in  this  problem  is  the  motion  of  a  surface  ship  at  sea 
during  an  aircraft  recover}’  phase.  This  is  a  linear,  time-angle  system  that  can  be  de¬ 
scribed  in  a  discrete  equation  of  ship  motion  .  The  state  space  equation  is 


**+i  =  ®k*k  +  r>* 


(2.1) 


where 

«£*  =  parameter  to  be  estimated  (state  vector)  heave,  pitch  and  roll 

<Ih=  state  transition  matrix  which  describes  how  the  states  of  the  dynamic  system  are 
related. 

T*=  system  random  uncorrelated  coefficient  matrix. 
uk-  sequence  of  random  uncorrelated  inputs. 

From  Equation  (2.1)  and  the  above  assumptions,  the  state  vector  for  heave,  pitch  and 
roll  is 


£k  = 


(2.2) 


and  the  system  state  equation  can  be  approximated  for  T  <  <  1 
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where 

uA  =  £„cosav  +  plant  noise 

C„  =  swell  wave  amplitude 

cos=  swell  angular  frequency 

c2  =  g,=  heaving  restoring  force  coefficient 

b2  =  heaving  damping  force  coefficient. 

cp  =  gf~  pitching  restoring  moment  coefficient 

bf  =  pitching  damping  moment  coefficient 

c,  =  g,  =  rolling  restoring  moment  coefficient 

br=  heaving  damping  moment  coefficient. 


The  system  noise  process  for  the  ship  motion  estimation  problem  is  a  function  of  the 
random  uncorrelated  coefficient  matrix.  Tk ,  and  the  random  forcing  function,  wA. 

D.  MEASUREMENT  MODEL 

For  a  linear  measurement  process,  the  measurements  are  linearly  related  to  the  state 
variables  and  can  be  modeled  using  the  linear  measurement  equation 


2-k  ~  +  h 


(2.4) 


where 

2A  =  set  of  measurements. 

II,  -  observation  matrix  that  gives  the  noiseless  relationship  between  the  measure¬ 
ments  and  the  state  vector. 

£,=  state  vector 

v,-  sequence  of  random  uncorrelated  measurement  noise. 

In  this  ship  motion  estimation  problem,  the  measurements  are  the  heave,  pitch  and 
roll  angles  taken  from  the  ship  sc  sors. 
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Ik  =  Cl  o  +  vk  (2.5) 

The  linear  equation  with  the  measurement  data  available  can  now  be  processed  with 
the  Kalman  filter.  The  type  of  noise  that  affect  the  ship  motion  data  are  assumed  to  be 
white  gaussian  noise.  This  noise  model  was  used  in  the  computer  simulations  and  shown 
in  Figure  2,  with  zero  mean  and  variance  of  one. 

The  complete  development  of  this  model  can  be  found  in  Bhattacharyya  [Ref.  1: 
pp.35-lS0J,  Walt  [Ref.  2:  pp. 31-47],  Korvin-  Kroukovsky  [Ref.  3]  and  Blagoveshchensky 
[Ref.  4], 

The  corresponding  block  diagram  of  the  system  is  shown  in  Figure  3. 


Figure  2.  White  Noise  Model 
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Figure  3.  Ship  Motion  Model  Block  Diagram 


III.  KALMAN  FILTER  THEORY 

A.  GENERAL 

Filtering  refers  to  the  process  of  estimating  the  state  vector  at  the  current  time  based 
upon  all  past  measurements.  An  optimal  filter  concentrates  on  optimizing  a  specific 
performance  measure  used  to  determine  the  quality  of  the  estimate.  The  Kalman  filter 
is  an  optimal  filter  in  a  class  of  linear  filters  that  minimize  the  mean  square  estimation 
error  between  the  actual  and  desired  output.  In  other  words,  the  Kalman  filter  attempts 
to  minimize  the  elements  along  the  main  diagonal  of  the  state  error  covariance  matrix. 
The  filter  is  a  recursive  algorithm  for  processing  discrete  measurements  or  observations 
in  an  optimal  manner,  [Ref.  5:  p.  101)  and  [Ref.  6).  It  requires  a  priori  knowledge  of  the 
state  estimate  (£*_,)  and  its  error  covariance  (/\_,),  and  the  current  observation  (ik).  The 
Kalman  filter  is  the  proper  algorithm  to  be  used  when  both  the  system  model  and  the 
measurement  model  are  linear  functions  of  the  state  variables  and  these  models  can  be 
described  by  the  equations 


£k+l  =®k£k  +  rkUk 

(3.1) 

±k  ~  ^hr—k  4"  ~k 

(3.2) 

In  this  equation,  ®k  and  T*  are  constant  matrices  and  H*  is  a  linear  constant  func¬ 
tion  of  the  state  variable  while  is  the  random  forcing  function  and  i*  is  the  meas¬ 
urement  noise.  1  he  plant  random  forcing  function  and  measurement  noise  are  assumed 
to  be  uncorrclated  (white  Gaussian)  with  zero  mean.  That  is 

F[vv(/c).wr(/)]  =  Q'(k)Skj  (3.3) 

and 


£[v(£).v7(/)]  =  R'{k)5kJ 


(3.4) 


where 


\  =  1,  k  —j 
=  0.  k  j*  j 


(3.5) 
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B.  NOISE  PROCESSES 

The  calculation  of  the  error  covariance  matrix  and  the  filter  gain  matrix  requires  the 
covariance  matrices  for  the  uncorrelated  noise  process  uk  and  For  the  measurement 
noise  process,  v*,  the  covariance  matrix  is 

-  p~k  (3-6) 

where  Rk  is  defined  as  the  state  measurement  noise  covariance  matrix.  It  is  based  on  the 
sensor  accuracy  and  accounts  for  unknown  disturbances  such  as  steps,  white  noise,  or 
imperfections  in  the  plant  model.  The  variance  of  the  white  noise  model  used  in  the 
computer  simulations  was  one. 

The  state  excitation  matrix,  Qk,  used  in  the  Kalman  filter  represents  the  system  noise 
process  and  is  a  function  of  the  system  noise  coefficient  matrix,  T*,  and  the  random 
forcing  funcaon,  uk.  This  matrix  is  given  by 

Qk  =  U'kQ‘krTk]  (3.7) 

where  T*  is  the  same  as  in  Equation  (2.3).  The  Qk  matrix  allows  for  any  random  heave, 
pitch  and  roll  ship  motion  well  as  inaccuracies  in  the  system  model.  The  magnitude  of 
Ok  has  a  direct  bearing  on  the  magnitude  of  the  state  error  covariance  matrix  and  it 
prevents  the  covariance  matrix  from  becoming  singular  by  ensuring  some  uncertainty  in 
the  state  estimates. 

C.  INITIALIZATION  AND  OPERATION 

In  the  ship  motion  estimation,  Kalman  filters  are  used  to  minimize  the  ship's  atti¬ 
tude  errors.  Prior  to  processing  the  measurement  data,  the  filter  must  be  initialized  with 
an  initial  state  estimate  and  an  initial  error  covariance  matrix.  This  initialization  process 
is  a  very  important  step  in  the  filter  operation  and  gross  inaccuracies  in  this  step  may 
cause  the  filter  to  diverge.  Divergence  occurs  when  the  calculated  covariance  errors  be¬ 
come  much  smaller  than  the  actual  covariance  errors.  This  causes  the  estimate  values 
of  the  states  to  pull  away  from  the  actual  value  and  is  also  used  as  a  startup  value  for 
the  recursive  scheme. 

The  basic  operation  of  the  filter  is  a  relatively  straightforward  recursive  process. 
Based  on  Kirk  (Ref.  6  :  pp.4-1,  4-62]  and  Gelb  [Ref.  5J.  the  discrete  time  version  of 
equations  used  in  the  Kalman  filter  are: 

£(k\k-\)  ~  ^k^{k\k)  (3.8) 
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R(k\k- 1)  =  ^k^iklk^k  +  Qk 

(3.9) 

Gk  =  P(k\k-\)Hk(HkR(k\k-\)Hk  +  &k) 

(3.10) 

&{k\k)  =  &{k\k-\)  +  Gk(zk  -  (*_!)) 

(3.11) 

R(k\k)  =  (I  ~  GkHk)P(k\k- 1) 

(3.12) 

where 

£**_,>  =  projected  ahead  state  estimate 

4>t  =  state  transition  matrix  given  by  Eq.  (2.3) 

P(**_i>  -  projected  ahead  state  error  covariance  matrix 
Qk  =  state  excitation  covariance  matrix  given  by  Eq.  (3.4) 

Gk  =  Kalman  gain  matrix 

Rk  =  state  measurement  noise  covariance  matrix  given  by  Eq.  (3.3) 

Hk  =  linearized  measurement  matrix  given  by  Eq.  (2.5). 

From  the  given  equations,  we  can  start  the  filter  processing  operation.  The  a  priori 
state  estimate  are  calculated  using  the  O  matrix  shown  in  Equation  (2.3)  where  7' is  the 
sampling  time  in  seconds  between  the  observed  ship  motions  (heave,  pitch  and  roll). 
(It  is  assumed  ship  motions  arc  received  simultaneously  through  the  ship  sensors.) 

The  Kalman  gain  matrix  serves  to  minimize  the  mean  square  estimation  error  and 
is  an  indication  of  how  much  emphasis  or  weight  will  be  placed  on  the  current  observa¬ 
tion.  If  is  small,  the  Kalman  gain  matrix  will  also  be  small  due  to  the  finite  value 
of  Rk.  If  the  P{kk. „  is  relatively  large,  the  gain  is  approximately  one.  By  rewriting  the 
equation  for  the  calculation  of  the  state  estimate  Equation  (3.11).  as 

£{k\k)  =  (I  ~  GkHk)X{k\k-\)  +  %|(  (3.13) 

we  can  see  how  the  Kalman  gain  matrix  directly  affects  the  weight  placed  on  the  current 
observation  zk.  A  large  gain,  indicating  a  large  a  priori  error  covariance,  will  place  more 
weight  on  the  current  observation  as  the  filter  tries  to  correct  the  states.  A  small  gain, 
indicating  a  small  error  covariance,  places  less  emphasis  on  the  new  observation. 

If  the  Kalman  gain  is  expressed  as 

~  P{k\k~])^k‘ Rk  '  (3.14) 

i 
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it  can  be  seen  that  the  gain  matrix  is  "proportional"  to  the  uncertainty  in  the  estimate 
and  "inversely  proportional"  to  the  measurement  noise  Rk.  For  a  large  R„  and  a 
small  /’(A*-!,,  the  measurement  in  Equation  (3.2)  is  due  mainly  to  noise  and  only  small 
corrections  should  be  made  in  the  state  estimate.  However,  if  Rk  is  small  and  P(4*_n  is 
large,  the  measurement  contains  considerable  information  about  the  errors  in  the  esti¬ 
mates  and  therefore,  a  strong  correction  should  be  made  to  the  state  estimates  [Ref.  5: 
pp.  127-8], 

The  computation  process  for  prediction  is  divided  into  the  following  steps:  1)  cal¬ 
culate  the  estimate  using  Equation  (3.11)  and  2)  use  Equation  (3.8)  for  the  desired  pre¬ 
diction. 

The  block  diagram  of  the  system  with  Kalman  Filter  is  shown  in  Figure  4  and  the 
system  predictor  block  diagram  is  shown  in  Figure  5. 
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Figure  A.  Ship  Motion  Model  Diagram  With  Estimator 
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Figure  5.  Ship's  Motion  Predictor  Block  Diagram 


IV.  COMPUTER  SIMULATIONS 


A.  GENERAL 

A  measured  simulated  ship  motion  (heave,  pitch  and  roll)  and  the  forcing  function 
data  were  chosen  based  from  the  ship  speed  of  25  knots,  heading  000  degrees  true  and 
the  swell  wave  heading  direction  at  120  degrees  true.  The  sea  state  for  a  swell  wave 
height  of  11  feet  is  between  5  and  6  based  on  Bhattacharyya  [Ref.  1:  p.  104]  which  can 
be  described  that  large  waves  begins  to  form  while  white  crests  are  extensive  everywhere 
(probably  some  spray).  The  wind  force  (beaufort)  is  classified  as  level  6  which  is  char¬ 
acterized  by  a  strong  breeze  with  the  velocity  of  25  knots. 

These  measured  data  were  used  in  the  Matlab  computer  program  ([Ref.  7]  and  [Ref. 
8])  algorithm  (as  shown  in  Appendix  A)  to  generate  the  graphical  results  to  see  the 
physical  representation  of  the  ship  motion  oscillations.  These  graphical  representations 
of  swell,  heave,  pitch  and  roll,  which  are  the  results  from  the  above  computer  runs,  have 
the  following  parameters  listed  in  Table  1  and  are  shown  in  Figures  6,  7,  S  and  9. 


Table  1.  WAVE  AND  SHIP  MOTION  PARAMETERS,  TS  =  0.1  SEC 


Coefficient 

Swell 

Heave 

Pitch 

Roll 

Wave  ht  (ft) 

11 

— 

— 

-- 

Period  (sec) 

7 

11 

10 

18 

«y  (rad  sec) 

0.8976 

0.5712 

0.6283 

0.3491 

Damping 

0 

0.6S54 

0.7540 

0.4189 

Restoring 

0 

0.62S3 

0.3491 

Kalman  Gain 

- 

wmmmsm 

mmm 

0.1191 

0.0756 

mm 
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Figure  6.  Swell  Wave  Motion  with  a  Height  of  11  feet 
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Figure  8.  Ship's  Pitch  Motion  caused  by  11  feet  high  of  Snell 
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Figure  9.  Ship's  Roll  Motion  caused  by  11  feet  high  of  Snell 
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With  the  simulated  observed  ship  motion  data  for  heave,  pitch  and  roll  generated 
from  the  above  computer  runs,  the  Kalman  filter  was  used  to  calculate  the  estimates  by 
using  the  Matlab  computer  program  algorithm  generated  as  shown  in  Appendix  B.  The 
resulting  estimated  data  were  represented  in  graphical  form  to  show  the  difference  be¬ 
tween  the  observed  and  the  estimated  ship  motions.  These  graphs,  with  estimates,  are 
shown  in  Figures  10,  11  and  12.  Computer  runs  were  also  conducted  to  show  the 
graphical  representations  of  the  predicted  ship  motions  heave,  pitch  and  roll.  A  Matlab 
computer  program  is  shown  in  Appendix  C  for  this  calculation.  These  graphical  pred¬ 
ictions  of  ship  motions  are  shown  in  Figures  13,  14  and  15.  As  we  can  see,  the  estimates 
and  predictions  illustrate  a  well-defined  representation  of  the  ship  motion  that  can  be 
used  in  the  prediction  of  the  ship  motion  in  the  future  time  w'hich  is  to  be  relayed  to  the 
Landing  Signal  Officer  for  his  decision  for  waveoffs  and  to  the  aircraft's  pilot  through 
the  aircraft's  TACAN  during  the  recovery  phase  of  flight  operations. 


ps  Observed  Roll  and  Roll  Estimate 


Time(sec) 


Figure  15.  Observed  and  Predicted  Roll  Motion  caused  by  1 1  feet  high  of  Swell 
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V.  CONCLUSIONS 


The  feasibility  of  estimating  an  aircraft  carrier  motion  at  sea  by  measuring  the  ship's 
actual  position  and  motion  was  investigated.  The  ship's  motion  representation  was 
generated  based  on  ship's  motion  mathematical  model  taken  from  [Ref.  lj.  Subse¬ 
quently,  a  Kalman  filter-estimator  adapted  for  real-time  computation  on  a  digital  com¬ 
puter,  was  generated.  The  results  obtained  show  that  the  desired  prediction  time  can 
be  reached  with  reasonably  acceptable  errors. 

Being  able  to  predict  the  ship's  motion  can  lead  to  an  improvement  of  aircraft 
landing  accuracy  and  safety.  This  can  be  accomplished,  for  instance,  by  generating  ter¬ 
minal  guidance  (landing)  laws  making  use  of  the  future  ship's  motion  and  position. 
Moreover,  the  possibility  of  prediction  can  improve  the  LSO  information  and  policy  for 
landing  acceptance  or  wave-oils.  The  possibility  of  processing  the  measured  motion  by 
Fast  Fourier  Transform  Algorithms  (FFT),  in  order  to  obtain  the  estimated  ship  motion 
parameters  such  as  the  heave  (  z  ),  pitch  (8)  and  roll  (<£),  in  real  time  may  lead  toward 
an  adaptive  predictor  scheme. 
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APPENDIX  A.  MATLAB  PROGRAM  1 


This  Matlab  computer  program  calculates  the  Observed  Ship's  Motion  such  as  the 
heave,  pitch  and  roll  of  the  ship.  Required  coefficients  must  be  available  to  run  the 
program. 


%  A  Matlab  program  to  generate  the  ship's  heave,  pitch  and  roll 
X  roll  motion  with  a  given  Swell  wave  force  of  N  sin  ws*t. 

% 

% 

mmmmmmmfflH  swell  parameter  xxxxxxxxxxxxxxxxxxxxxxxxxxx 
% 

% 

% 

diary  p. m 

%  Swell  period  (Ts)  seconds 
Ts=7; 


% 


Swell  angular  frequency 

ws=  2*pi*l/Ts  ; 

Wave  Height  (N)  feet 
N-ll; 
bl=N*ws; 


% 

X 


Restoring  coefficient 
cs=ws**2 


Heave  parameter  XXXXXXXXXXXXXXXXXXXXXXXXXXX 


0/  Of  Of  O/  Of  Of  Of  Of  0/  Of  0/  Of  Of  Of  Of  Of  Of  Of  Of  Of  0/  0/  Of  0/  o.' 

/o  /o  /o  /o  /o  -  o  fo  /o  /o  /o  /o  /g  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  o 

% 

%  Heave  period  (th)  seconds  typical  for  carriers 
th=  11; 


X 

X 


Heave  angular  frequency  (wh)  rad/sec 
wh=  2*pi*l/th  ; 


Heave  restoring  cefficient  (cz) 
cz=wh 


7o  Heave  damping  factor  (zetah) 
zetah=  0. 6; 


%  Damping  coefficient  (bh) 
bz=2*zet  ah*wh 

Of 

/o 


X 

Of  Of  0/  O/  0/  Of  0/  O/  0/0/ Of  Of  0/  O/  0/  0/  0/  0/  0/  Of  0/  0/  O,'  O/  0/ 

/o  'O  /o  /o  /o  /o  /o  /O  lo'olofo  >0  to  fo  to  /O  /O  /O  /O  /o  /O  'O  /o  o 

O' 

fo 


Pitch  parameter 


Of  Of  Of  0/  O'  Of  Of  Of  Of  Of  0>  Of  Of  Of  0/  Of  Of  Of  Of  Of  Of  Of  Of  Of  0  ■  0/  Of 

fo  fo  fo  fo  fo  fo  fo  'O  fo  fo  >0  fa  'O  fa  fa  fa  fo  fo  fa  >0  fa  fa  >o  fa  to  fa  fa 


[ 


s>?  s?  o-?  a-s  o'?  a-?  a-?  a-?  a-?  a^a^ae  a?  a*a« 


Pitch  period  (tp)  seconds  typical  for  carriers 
tp=10; 


Pitch  angular  frequency  (wp)  rad/sec 
wp=  2*pi*l/tp  ; 


Pitch  Restoring  coefficient 
cp=wp 

Heave  damping  factor  (zetap) 
zetap=  0.  6 

Damping  coefficient  (bp) 
bp=2*zetap*wp 


Roll  parameter  XX%%XXXX%%XXX%%%%%\%%%%%%%%%% 


Roll  period  (tr)  seconds  typical  for  carriers 
tr=  18; 

Roll  angular  frequency  (wr)  rad/sec 
wr=  2*pi*l/tr  ; 


Roll  Restoring  coefficient 
cr=wr 

Roll  damping  factor  (zetar) 
zetar=  0. 6  ; 

Damping  coefficient  (br) 
br=2,vzetar,vwr 


if  0/0/ Of  Of  Of  O/  Of  O'  0/  O/  O/  O/  0/  0/  ©/  Of  0/  Of  0 /  Of  0/ 

o/o/o  /o/o/o/o/o  c  to  /©  >o  Jo  /o  /o  /o  >o  to  to  /o 


Parameters  of  the  system 


where: 

a=[];--  the  A  martix  coefficient  in  continous  time 
b=[ ]  ;  -  -  the  B  matrix  coefficient  in  continuos  time 


Swell  Parameters 


as=-0  1; -cs  0-; 
bs=[  0;  1]  ; 


Heave  Parameters 


ah=[0  1;  -cz**2  -bz]  ; 
bh=[  0; .  3] ; 

% 

%  Pitch  Parameters 

%  ======= 

% 

ap=[ 0  1; -cp**2  -bp]  ; 
bp=[  0; .  3] ; 

% 

% 

7o  Roll  Parameters 

%  ============== 

% 

ar=[  0  1;  -cr*,v2  -br]  ; 
br=[ 0; . 6] ; 
diary  off 
% 


%  Measurement  Parameters 

%  ======= == 

7<  c=[  1  0]; 

d=[  0]  ; 
dt=0.  1; 
q=eye( 2); 
g=eye(2); 
r=l/10000; 
h=[  1  0] ; 
cs=[ N  0] ; 

% 

%%%%%%%%%%%%  Determine  the  discrete  coeficient  of  the  system  %%%%%%%%% 
%  The  Phi's  and  Gamma’s  of  the  system 

Of 

/O 

0/ 

/o 

[ phih.gamh] =c2d(ah,bh ,dt)  ; 

[phip.gamp] =c2d(ap,bp,dt)  ; 

[phir,gamr]=c2d(ar,br,dt)  ; 

[ phis , gams] =c2d( as ,bs ,dt)  ; 

% 


%%*%%%%%%%%%  Generate  Input  for  the  system  %%%%%%%%%%%%%%%%%%%%%%%%%% 


t=( 0: . 1: 100)'; 
rand( 'normal' ); 
v=sqrt(r)*rand( 1000, 1); 


us=dlsim(phis ,gams ,cs ,d,v); 


%  Time  axis 

%  Random  are  uniformly  distributed 
%  in  the  interval  (0. 0,1.0).  Random 
%  ('normal')  switches  to  normal  with 
7e  mean  of  0  and  variance  of  1 


%%%%%%%%%%%%  Generate  output  values  for  the  system  %%%%%%%%%%%%%%%%%%% 

?. 


yh=d Is im( phih , gamh , c , d , us ) ; 
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yp=d Is im( phip , gamp , c , d , us ) ; 
yr=dlsim(phir ,gamr,c,d,us); 


%%%%%%%%%%%%%  Plot  the  generated  output  vs  time  of  the  system  %%%%%%% 

% 

plot(yh) ,title( ' Ships  Heave  Motion'); 
xlabel('Time  (sec) ’) ,ylabel( 'Heave  angle  (rad)'); 
delete  pitch. met 
meta  pitch 

plot(yp) ,title( ' Ships  Pitch  Motion'); 

xlabel( 'Time  (sec) ') ,ylabel( 1  Pitch  angle  (rad)'); 

meta 

plot(yr) ,title( ' Ships  Roll  Motion'); 

xlabel( 'Time  (sec) ') ,ylabel( 'Roll  angle  (rad)'); 

meta 

plot(v) .title( 'Random  Noise  vs  Time'); 

xlabel( 'Time  (sec) ') ,ylabel( ' Noise  Magnitude'); 

meta 

plot(us) ,title( ' Swell  with  White  Noise  vs  Time'); 
xlabel( 'Time  ( sec) ’) ,ylabel( ' Swell  angle  (rad)'); 
meta 
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APPENDIX  B.  MATLAB  PROGRAM  2 


This  Matlab  computer  program  calculates  the  Observed  Ship's  Motion  and  esti¬ 
mates  in  real  time  such  as  the  heave,  pitch  and  roll  of  the  ship.  Required  coefficients 
must  be  available  to  run  the  program. 

%  This  Matlab  computer  program  calculates  the  Observed  and  Estimated 
%  state  of  the  Ship  s  Heave,  Pitch  and  Roll  based  on  the  given  data. 

% 


%%%% %%%%%%%%%%%%%%%%%%%%%%  SWELL  parameter 

% 


0/0/ 0/0/ Of  0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 
/oloioiototo/otohloio/ohfololotolotolotolotoloiololo 


diary  p. m 

%  Swell  period  (Ts)  seconds 
Ts=7  ; 


%  Swell  angular  frequency 

ws=  2*pi*l/Ts  ; 


%  Wave  Height  (N)  feet 
N=ll; 


%  Restoring  coefficient 
cs=ws**2  ; 

% 

% 


Heave  parameter 


OJ  0/0/ 0/0/  O/  O/  O/  O/  O/  0/  O/  0/  0/  0/  0/  O/  Of  0/  O/  Of  O/  Oj  0/  0/  o/  o/ 
/o/o/o  fo  /o  /o  /o  /o  /o  /o  /o  / 6  /o  to  /o  /o  /o  /o  /o  /o  /o  fo  /o  /o  /o  to  to 


O' 

!o 


Heave  period  (th)  seconds  typical  for  carriers 
th=  11  • 

Heave  angular  frequency  (wh)  rad/sec 
wh=  2*pi*l/th  ; 


% 

%  Heave  restoring  cefficient  (cz) 
cz=wh  ; 

%  Heave  damping  factor  (zetah) 
zetah=  0. 6; 

%  Damping  coefficient  (bh) 

bz=2*zetah*wh  ; 


% 

% 

Pitch  parameter  %%%%%%%%%%%%%%%%%%% 

% 

% 

%  Pitch  period  (tp)  seconds  typical  for  carriers 


0/  O/  0/  Of  O/OJ  0/0/ 

/o  fo  fo /o/o/o  to  fo 


tp=10; 


%  Pitch  angular  frequency  (wp)  rad/sec 
wp=  2*pi*l/tp  ; 

% 

% 

%  Pitch  Restoring  coefficient 

cp=wp  ; 

% 

%  Heave  damping  factor  (zetap) 

zetap=  0.  6  ; 

% 

%  Damping  coefficient  (bp) 

bp=2*zetap*wp  ; 


Roi  1  parameter  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 

%  Roll  period  (tr)  seconds  typical  for  carriers 
tr=  18; 


%  Roll  angular  frequency  (wr)  rad/sec 
wr=  2*pi*l/tr  ; 


%  Roll  Restoring  coefficient 

cr=wr 

% 

%  Roll  damping  factor  (zetar) 

zetar=  0. 6  ; 

%  Damping  coefficient  (br) 

br=2*zetar*wr  ; 

•/ 

'O 

O/ 

/o 

0/  O/  0/  O/  O/  O/  0/  0/  O/  0/  O/  Of  0/  0/  0/  Of  Of  O/  0/  O/  Of  O/  0/  O/  Oy 


%%%%%%%%%%%%%%%%% %%%%%%%%  Parameters  of  the  system  %%%%%%%%%%%%%%%%%%%% 
% 


where: 


/0 

% 

0/ 

/o 

% 


a=[ ] ; --  the  A  martix  coefficient  in  continous  time 
b=[ ] ; —  the  B  matrix  coefficient  in  continuos  time 

Swell  Parameters 


as=-0  1; -cs  0-; 
fs=[  0;  1] ; 


Heave  Parameters 


ah=[  0  1; -cz**2  -bz]  ; 
fh=[ 0; .  3]  ; 
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sa*s$3'?  a«  a?  a-s  as  as  a-s  a^s  as  as  as  as  as a? as  as  as  as  ^ i? ^ ^ s< a5  q.  as  a*  as  as  as  as  as  as 


Pitch  Parameters 


ap=[ 0  1; -cp**2  -bp]  ; 
fp=[  0; .  3]  j 


Roll  Parameters 


ar=[ 0  1; -cr**2  -br]  ; 

fr=[ 0; -  6] ; 

Lary  off 


Measurement  Parameters 


c=[  1  0]; 
dt=0.  1; 

Swell  Noise 

q=l/100; 

Measurement  Noise 
r=l/8500; 


d=[  0]  ; 

Determine  the  Discrete  coeficient  of  the  system 
Swell 

[  phis  , gams]  =c2d(as  ,  fs  ,dt)  ; 

HEAVE 

[phih,gamh]=c2d(ah,fh,dt)  ; 

PITCH 

[ phip ,gamp] =c2d( ap , fp ,dt)  ; 

ROLL 

/O 

[  phir ,gamr] =c2d(ar , fr ,dt)  ; 

% 


0/0/  0/0/ 
/o/o>o/o 
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%XX%\%%%X%  Calculate  the  Kalman  Gains  for  Heave,  Pitch  and  Roll  XXXXXXX 
X 

kh=dlqe(phih,gamh,c,q,r)  ; 
kp=d lqe( ph ip , gamp ,  c ,  q ,  r  )  ; 
kr=dlqe(phir ,gamr,c,q,r)  ; 

% 

XXXXXXXXXXXXXXXXXXXXXX  Initial  Condition  XXXXXXXXXXXXXXXXXXXXXXX 
% 

xh=[  0.  05;  0.  05]  ; 
xhe=[  0.05;  0.05]  ; 

xp=[  0.  05;  0.  05]  ; 
xpe=[ 0. 05; 0. 05]  ; 

xr=[0.  050;  0.0506]  ; 
xre=[  0.  05; 0. 05]  ; 
xs=[  0.  05;  0.  05]  ; 

% 

XXXXXXXXXXXX XXXXXXXXX  Random  Noise  Generators  XXXXX XXXXXXXXXXXXX 


rand( 'normal' ) 
w=sqrt(q)*rand(  1 , 1500); 
v=sqrt( r)*rand( 1 , 1500); 
t=0:  .  1:  100; 

% 

XXXXXXXXXXXXXXXXXXXX  Simulation  Loop  XXXXXXXXXXXXXXXXXXXXXXXXX 
% 

for  k=l:  1000; 

% 

XXXXXXXXXXXXXXXXXXXX  Input  for  Heave,  Pitch  nad  Roll  %%%%%%%%% 

o / 
to 

xs(: ,k+l)=phis*xs( :  ,k)  +  gams*w(k); 

X 

XXXXXXXXX  State  calculation  for  Heave,  Pitch  and  Roll  "XXXXXXXXX 

O/ 

to 

xh( : ,k+l)=phih*xh( :  ,k)  +  gamh*xs(k)  ; 
xp(:  ,k+l)=phip*xp( :  ,k)  +  gamp*xs(k)  ; 
xr(: ,k+l)=phir*xr( :  ,k)  +  gamr*xs(k)  ; 


XXXXXXXXX  Output  calculation  for  Heave,  Pitch  and  Roll  in  Discrete  form 

O' 

yh(:  ,k+l)=c*xh( :  ,k+l)  +  v(k+l); 
yp( :  ,k+l)=c*xp(:  ,k+l)  +  v(k+l); 
yr(:  ,k+l)=c*xr(:  ,k+l)  +  v(k+l); 


XXXXXXXXX  State  estimate  calculation  for  Heave,  Pitch  and  Roll  XXXXXXXXXX 
7 

xhe(: ,k+l)=  phih*xhe(:  ,k)  +  kh*(yh(:  ,k+l) -c*(phih*. . . 
xhe(: ,k)  )); 

xpe( :  ,k+l)=  phip*xpe(:  ,k)  +  kp*(yp(:  ,k+l)-c*(phip*.  .  . 
xpe(: ,k)  )); 

xre(:  ,k+l)=  phir*xre(:  ,k)  +  kr*(yr(:  ,k+l) -c*( phir*.  .  . 
xre(: ,k)  )); 
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end 


plot(t,xh( 1,:  ) ,t,xhe( 1 ))  ,title( ' Ships 
Observed  Heave  and  Heave  Estimate') 

xlabel( 'Time(sec) ') ,ylabel( 'Heave  Angle  (rad)') 
delete  esth.  met 
meta  est 

plot(t,xp( 1,: ) ,t ,xpe( 1 , : )) ,title( ' Ships 
Observed  Pitch  and  Pitch  Estimate') 

xlabel( 'Time(sec) ') ,ylabel( ' Pitch  Angle  (rad)') 
meta 

plot(t ,xr( 1,: ) ,t ,xre( 1,:  )) ,title( ' Ships 
Observed  Roll  and  Roll  Estimate') 

xlabel( 'Time(sec) ' ) ,ylabel( 'Roll  Angle  (rad)') 
meta 
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APPENDIX  C.  MATLAB  PROGRAM  3 

This  Matlab  computer  program  calculates  the  Observed  Ship's  Motion  and  predicts 
in  real  time  such  as  the  heave,  pitch  and  roll  of  the  ship.  Required  coefficients  must  be 
available  to  run  the  program. 


%  This  Matlab  computer  program  calculates  the  Observed  and  Predicted 
%  state  of  the  Ship  s  Heave,  Pitch  and  Roll  based  from  the  given  data. 
% 


% 

% 

SWELL  parameter 

1 

% 

o/ 

/o 

%  Swell  period  (Ts)  seconds 
Ts=7  ; 


7o  Swell  angular  frequency 

ws=  2*pi*l/Ts  ; 


%  Wave  Height  (N)  feet 

N=ll; 

%  Restoring  coefficient 
cs=ws**2  ; 

% 

o/ 

ramsmHKKnaa  Heave  parameter 

% 

%  Heave  period  (th)  seconds  typical  for  carriers 

th=  11  ; 

O/ 

A> 

%  Heave  angular  frequency  (wh)  rad/sec 

wh=  2*pi,vl/th  ; 


% 

%  Heave  restoring  cefficient  (cz) 

cz=wh 


%  Heave  damping  factor  (zetah) 
zetah=  0. 6; 


%  Damping  coefficient  (bh) 

bz=2*zetah,vwh  ; 

% 

7o 

Pitch  parameter  %%%%%%%%%%%%%%%%%%%%%%%%%%% 

0/ 


%  Pitch  period  (tp)  seconds  typical  for  carriers 
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tp=10; 


%  Pitch  angular  frequency  (wp)  rad/sec 

wp=  2*pi*l/tp  ; 

% 

% 

%  Pitch  Restoring  coefficient 

cp=wp  ; 

% 

%  Heave  damping  factor  (zetap) 

zetap=  0. 6  ; 

% 

%  Damping  coefficient  (bp) 

bp=2*zetap*wp  ; 

% 

% 

mmmmnmmma  Ron  parameter  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 

%  Roll  period  (tr)  seconds  typical  for  carriers 

tr=  18; 

%  Roll  angular  frequency  (wr)  rad/sec 

wr=  2*pi*l/tr  ; 


%  Roll  Restoring  coefficient 

cr=wr 

% 

%  Roll  damping  factor  (zetar) 
zetar=  0. 6  ; 


Damping  coefficient  (br) 

br=2*zetar*wr 


% 

% 

o/  0/  O/  0/  O/  0/  O/  0/  0/  0/  0/  0/  0/  0/  0/  O/  0/  O/  O/  O/  0/  O/  O/  0/  0/ 
/o  /o  /o  /o  /o  /Q  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /©  /o  /o  /o  /o  /o  /o 

% 


Parameters  of  the  system  %%%%%%%%%%%%?£%%%%%%% 


%  where: 

%  a=[ ] ; —  the  A  martix  coefficient  in  continous  time 

%  b=[];--  the  B  matrix  coefficient  in  continuos  time 

% 


% 

0/ 

/o 

% 


% 

% 

% 


Swell  Parameters 


as=[  0  1;  -cs  0]  ; 
fs=[  0;  1]  ; 


Heave  Parameters 


% 

ah=[ 0  1; -cz**2  -bz]  ; 
fh=[ 0; .  3]  ; 

% 
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a?  a?  a?  a-?  a?  a?  a^  a?  as  a?  a?  a?  a?  a?  a?  a?  a?  a?  a?  a?  a?  a?  a?  a?  S'?  a?  a?  a?  a?  j?  a?  a?  as  as  a?  a?  as  a?  as 


Pitch  Parameters 


ap=[ 0  1; -cp**2  -bp]  ; 
fp=t  0; .  3] ; 

Roll  Parameters 


ar=[ 0  1; -cr**2  -br]  ; 

f r=[ 0; .  6] ; 


Measurement  Parameters 


c=[  1  0]; 

dt=0.  1; 

Swell  Noise 

q=l/100; 

Measurement  Noise 
r=l/3500; 

d=[  0]  ; 


Uo%%*V/0%%%%%%%%%  Determine  the  Discrete  coeficient  of  the  system  %%%% 

Swell 

[phis, gams] =c2d( as, fs,dt)  ; 

HEAVE 

[ phih.gamh] =c2d(ah, fh,dt)  ; 

PITCH 

[  phip , gamp] =c2d( ap , f p , dt )  ; 

ROLL 

[phir,gamr]=c2d(ar,fr,dt)  ; 

% 
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% 

%%%%%%%%%%  Calculate  the  Kalman  Gains  for  Heave,  Pitch  and  Roll  %%%%%%% 

% 

kh=dlqe(phih,gamh,c,q,r)  ; 
kp=d lqe( phip , gamp ,  c ,  q ,  r  )  ; 
kr=dlqe(phir,gamr ,c,q,r)  ; 

% 


nummmmmm  Initial  Condition  %%%%%%%%%%%%%%%%%%%%%%% 

% 

xh=[ “0.  05; -0. 05]  ; 
xhe=[ -0. 05; -0. 05]  ; 

xp=[  0.  05;  -0.  05]  ; 
xpe=[ -0.  05; -0. 05]  ; 

xr=[0.  050; -0.  050]  ; 
xre=[ -0. 05; -0. 05]  ; 

xs=[0.  05; -0.  05]  ;  ; 

% 

%%%%%%%%%%%%%%%%%%%%%  Random  Noise  Generators  XXXXXXXXXXXXXXXXXX 
% 

rand( 'normal' ) 
w=sqrt(q)*rand( 1,1200); 
v=sqrt(r)*rand( 1,1200); 
t=0:  .  1:  95; 

% 


0/  0/  0/  Of  Of  Of  0 f  Of  Of  Of  Of  Of  Of  Of  Of  Of  Of  Of  Of  Of 
fo  /o  to  to  fo  fo  fo  fo  to  to  to  fo  to  fo  fo  fo  fo  /o  fo  fo 

% 


Simulation  Loop  To  fo  /o% /o /o  -o  %  -4/0  fo%%7o%%%  /o%%  lo/ofo 


for  k=l:  950; 

% 


%%%%%%%%%%%%%%%%%%%%  Input  for  Heave,  Pitch  nad  Roll  %%%%%%%%% 


% 

o/ 

/o 


xs(:  ,k+l)=phis*xs(:  ,k)  +  gams*w(k); 


%%7o%%%%%%  State  calculation  for  Heave,  Pitch  and  Roll  %%%%%%%%% 

Of 


xh(:  ,k+l)=phih*xh(:  ,k)  +  gamh*xs(k)  ; 
xp(:  ,k+l)=phip*xp( :  ,k)  +  gamp*xs(k)  ; 
xr(: ,k+l)=phir*xr( :  ,k)  +  gamr*xs(k)  ; 


% 

1a7o%7o%°a  Output  calculation  for  Heave,  Pitch  and  Roll  in  Discrete  form 

% 


yh(:  ,k+l)=c*xh(:  ,k+l)  +  v(k+l); 
yp(:  ,k+l)=c*xp(: ,k+l)  +  v(k+l); 
yr(:  ,k+l)=c*xr(: ,k+l)  +  v(k+l); 

% 

% 

%%%%%%%  State  estimate  calculation  for  Heave,  Pitch  and  Roll  %%%%%%%%%% 

% 

xhe(:  ,k+l)=  phih*xhe(:  ,k)  +  kh*(yh(:  ,k+l)-c*(phih*.  .  . 
xhe( :  ,k)  )); 

xpe(:  ,k+l)=  phip*xpe(:  ,k)  +  kp*(yp(:  ,k+l)-c*(phip*. . . 
xpe(: , k)  )); 

xre(:  ,k+l)=  phir*xre( :  ,k)  +  kr*(yr(: ,k+l)-c*(phir*.  .  . 
xre(: ,k)  )); 

% 

% 
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a? 


% 

XXXXXX  State  Prediction  Calculation  for  Heave,  Pitch  and  Roll  XXXXXXX 

% 

xhp(:  ,k+l)=phih*xhe(: ,k)  +  gamh*xs(k)  ; 
xpp(:  ,k+l)=phip*xpe(:  ,k)  +  gamp,vxs(k)  ; 
xrp(:  ,k+l)=phir*xre(: ,k)  +  gamr*xs(k)  ; 

end 


plot(t,xh(l,: ,t,xhp(l,:  )) ,title( 1  Ships 
Observed  Heave  and  Predicted  Heave') 

xlabel( 'Time(sec) ') ,ylabel( 'Heave  Angle  (rad)') 
delete  predicth. met 
meta  predicth 

plot(t,xp( 1,: ),'-' ,t,xpp(l,: )) ,title( 'Ships 
Observed  Pitch  and  Predicted  Pitch') 

xlabel( 'Time(sec) ') ,ylabel( 1  Pitch  Angle  (rad)') 
meta 

plot(t ,xr( 1 ) , ' - ' ,t ,xrp( 1 ) ) ,title( ' Ships 
Observed  Roll  and  Predicted  Roll') 

xlabel( 'Time(sec) ' ) ,ylabel( 'Roll  Angle  (rad)') 
meta 
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